%%  
% R   = zeros(EEG.nbchan, EEG.pnts);
% b0  = zeros(EEG.nbchan, EEG.pnts);
% b1  = zeros(EEG.nbchan, EEG.pnts);
% b2  = zeros(EEG.nbchan, EEG.pnts);
% GL  = [ [ones(1,43) zeros(1,51)]' [zeros(1,51) ones(1,43)]'  ones(94,1) ];
% pGL = pinv(GL);
% 
% % compute R at all time and channels
% % ----------------------------------
% for c = 1:EEG.nbchan
%     for p = 1:EEG.pnts
%         tmp       = squeeze(EEG.data(c,p,:));
%         params    = pGL*tmp;
%         b0(c,p)   = params(3);
%         b1(c,p)   = params(1);
%         b2(c,p)   = params(2);
%         m         = corrcoef(GL*params, tmp);
%         R(c,p)    = m(2); %sqrt(mean((GL*params-tmp).^2));
%     end;
% end;
% 
% % find max R for all channels
% % ---------------------------
% for p = 1:EEG.pnts
%     [tmp mxind] = max(R(:,p));
%     R(EEG.nbchan+1,p) = tmp;
%     b0(EEG.nbchan+1,p) = b0(mxind,p);
%     b1(EEG.nbchan+1,p) = b1(mxind,p);
%     b2(EEG.nbchan+1,p) = b2(mxind,p);
% end;

%% test baseline for each subject
for d = 25 %1:3:length(ALLEEG)
    EEG1 = ALLEEG( d );
    EEG2 = ALLEEG( d+1 );
    EEG1.data = eeg_getdatact(EEG1, 'verbose', 'off');
    EEG2.data = eeg_getdatact(EEG2, 'verbose', 'off');
    
    naccu = 3000;
    R     = zeros(EEG1.nbchan, naccu);
    GL  = [ [ones(1,EEG1.trials) zeros(1,EEG2.trials)]' [zeros(1,EEG1.trials) ones(1,EEG2.trials)]'  ones(EEG1.trials+EEG2.trials,1) ];
    pGL = pinv(GL);

    range1 = [-2000 -1010];
    EEG1 = pop_rmbase(EEG1, range1);
    EEG2 = pop_rmbase(EEG2, range1);
    range2 = [-1000 -10];

    [tmp t1] = min(abs(EEG1.times-range2(1)));
    [tmp t2] = min(abs(EEG1.times-range2(2)));
    ERPmean  = [ squeeze(mean(EEG1.data(:,t1:t2,:),2)) squeeze(mean(EEG2.data(:,t1:t2,:),2)) ];
    
    % compute regular bootstrap
    
    ERPmean2 = { squeeze(mean(EEG1.data(:,t1:t2,:),2)) squeeze(mean(EEG2.data(:,t1:t2,:),2)) };
    [ t df pvalboot] = statcond( ERPmean2, 'mode', 'bootstrap', 'naccu', 1000);
    
    %ERPmean2(:,1:43)  = ERPmean1(:,1:43)+2;
    %ERPmean2(:,44:94) = ERPmean1(:,44:94)-2;

    % compute R at all time and channels
    % ----------------------------------
    for c = 1:EEG1.nbchan
        for ind = 1:naccu+1
            if ind == naccu+1
                tmp = ERPmean(c,:);
            else
                tmp = shuffle(ERPmean(c,:));
            end;
            params    = pGL*tmp';
            m         = corrcoef(GL*params, tmp);
            R(c,ind)  = m(2); %sqrt(mean((GL*params-tmp).^2));
        end;
    end;

    % find max R for all channels
    % ---------------------------
    [tmp maxE] = max(R(:, naccu+1));
    newR = max(R, [], 1);

    % compute bootstrap GLM
    % ---------------------
    orival       = newR(end);
    [surrog idx] = sort(newR);
    [tmp mx]     = max( idx );        

    len = length(surrog);
    pvals = 1-(mx-0.5)/len;
    pvals = min(pvals, 1-pvals);
    pvals = 2*pvals;
    fprintf('%s, %s, p-values: %2.4f (bootstrap %2.4f)\n', ALLEEG(d).subject, STUDY.datasetinfo(d).group, pvals, min(pvalboot));
end;

%% plot best result for richard
figure;plot(ERPmean(maxE,:));
x1 = 1:EEG1.trials;
x2 = EEG1.trials+1:size(ERPmean,2);
m1 = mean(ERPmean(maxE,x1));
m2 = mean(ERPmean(maxE,x2));
n1 = sum(ERPmean(maxE,x1) > (m1+m2)/2);
n2 = sum(ERPmean(maxE,x2) > (m1+m2)/2);
hold on; 
plot(x1, m1*ones(1,length(x1)), 'k');
plot(x2, m2*ones(1,length(x2)), 'k');
